﻿#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <windows.h>

// 配置数据文件存储目录
#define DATA_DIR "D:\\曲线拟合\\data\\"

// 函数声明
void mkdir(const char* path);
void generate_dataset(const char* filename, int t, double a, double b, double c, double d);
void least_squares_fit(const char* filename, int t);
void gradient_descent_fit(const char* filename, int t);
void visualize(const char* filename);

// 创建目录函数（带错误检查）
void mkdir(const char* path) {
    char cmd[512];
    sprintf(cmd, "mkdir \"%s\"", path);
    int result = system(cmd);

    // 检查目录是否存在或创建成功
    DWORD attr = GetFileAttributesA(path);
    if (attr == INVALID_FILE_ATTRIBUTES) {
        fprintf(stderr, "Failed to create directory: %s (error code: %d)\n", path, result);
        exit(EXIT_FAILURE);
    }
    printf("Directory check successful: %s\n", path);
}

// 生成数据集函数
void generate_dataset(const char* filename, int t, double a, double b, double c, double d) {
    printf("\n=== Generating dataset ===\n");

    // 构建完整路径
    char full_path[256];
    sprintf(full_path, "%s%s", DATA_DIR, filename);
    printf("Full path: %s\n", full_path);

    // 创建目录（如果不存在）
    mkdir(DATA_DIR);

    FILE* file = fopen(full_path, "w");
    if (!file) {
        perror("Error opening file for writing");
        exit(EXIT_FAILURE);
    }

    // 生成数据点
    printf("Generating data points...\n");
    for (double x = -100.0; x <= 100.0; x += 1.0) {
        double y;
        switch (t) {
        case 1: y = a * x + b; break;
        case 2: y = a * x * x + b * x + c; break;
        case 3: y = a * x * x * x + b * x * x + c * x + d; break;
        default:
            fprintf(stderr, "Invalid function type! t must be 1, 2, or 3\n");
            fclose(file);
            exit(EXIT_FAILURE);
        }
        fprintf(file, "%.2f %.2f\n", x, y);
    }
    fclose(file);
    printf("Dataset generated successfully in %s\n", full_path);
}

// 最小二乘法拟合函数
void least_squares_fit(const char* filename, int t) {
    printf("\n=== Starting least squares fitting ===\n");

    // 构建完整路径
    char full_path[256];
    sprintf(full_path, "%s%s", DATA_DIR, filename);
    printf("Trying to open file: %s\n", full_path);

    FILE* file = fopen(full_path, "r");
    if (!file) {
        perror("Error opening file for reading");
        printf("Hint: Check if the file exists in the directory: %s\n", DATA_DIR);
        return;
    }
    printf("File opened successfully\n");

    // 动态分配内存存储数据点
    int max_points = 2001;
    double* x = (double*)malloc(max_points * sizeof(double));
    double* y = (double*)malloc(max_points * sizeof(double));
    if (!x || !y) {
        perror("Memory allocation failed");
        fclose(file);
        return;
    }
    printf("Memory allocated for data points\n");

    // 读取数据点
    int n = 0;
    printf("Reading data points...\n");
    while (fscanf(file, "%lf %lf", &x[n], &y[n]) != EOF && n < max_points) {
        n++;
    }
    fclose(file);
    printf("Read %d data points\n", n);

    if (n < 2) {
        printf("Insufficient data points for fitting (need at least 2 points)\n");
        free(x);
        free(y);
        return;
    }

    // 根据函数类型进行拟合
    if (t == 1) { // 一次函数拟合 y = ax + b
        printf("Performing linear fitting (y = ax + b)...\n");

        double sum_x = 0, sum_y = 0, sum_xx = 0, sum_xy = 0;
        for (int i = 0; i < n; i++) {
            sum_x += x[i];
            sum_y += y[i];
            sum_xx += x[i] * x[i];
            sum_xy += x[i] * y[i];
        }

        double denominator = n * sum_xx - sum_x * sum_x;
        if (fabs(denominator) < 1e-10) {
            printf("Error: Division by zero (data points might be collinear)\n");
            free(x);
            free(y);
            return;
        }

        double a = (n * sum_xy - sum_x * sum_y) / denominator;
        double b = (sum_y - a * sum_x) / n;

        printf("Fitted parameters: a = %.6f, b = %.6f\n", a, b);
        printf("Fitted equation: y = %.4f * x + %.4f\n", a, b);

        // 计算剩余平方和和R²
        double sse = 0, sst = 0, mean_y = sum_y / n;
        for (int i = 0; i < n; i++) {
            double y_pred = a * x[i] + b;
            sse += pow(y[i] - y_pred, 2);
            sst += pow(y[i] - mean_y, 2);
        }
        double r_squared = 1.0 - (sse / sst);

        printf("Sum of squared errors: %.4f\n", sse);
        printf("Coefficient of determination (R²): %.4f\n", r_squared);
    }
    else if (t == 2) { // 二次函数拟合 y = ax² + bx + c
        printf("Performing quadratic fitting (y = ax² + bx + c)...\n");

        double sum_x = 0, sum_y = 0, sum_xx = 0, sum_xy = 0;
        double sum_xxx = 0, sum_xxxx = 0, sum_xxy = 0;

        for (int i = 0; i < n; i++) {
            double xi = x[i], yi = y[i];
            sum_x += xi;
            sum_y += yi;
            sum_xx += xi * xi;
            sum_xy += xi * yi;
            sum_xxx += xi * xi * xi;
            sum_xxxx += xi * xi * xi * xi;
            sum_xxy += xi * xi * yi;
        }

        // 构建线性方程组
        double A[3][4] = {
            {n,      sum_x,  sum_xx,  sum_y},
            {sum_x,  sum_xx, sum_xxx, sum_xy},
            {sum_xx, sum_xxx, sum_xxxx, sum_xxy}
        };

        // 高斯消元法解方程组
        for (int k = 0; k < 3; k++) {
            double pivot = A[k][k];
            if (fabs(pivot) < 1e-10) {
                printf("Error: Singular matrix encountered\n");
                free(x);
                free(y);
                return;
            }

            for (int j = k; j < 4; j++) {
                A[k][j] /= pivot;
            }

            for (int i = 0; i < 3; i++) {
                if (i != k) {
                    double factor = A[i][k];
                    for (int j = k; j < 4; j++) {
                        A[i][j] -= factor * A[k][j];
                    }
                }
            }
        }

        double c = A[0][3];
        double b = A[1][3];
        double a = A[2][3];

        printf("Fitted parameters: a = %.6f, b = %.6f, c = %.6f\n", a, b, c);
        printf("Fitted equation: y = %.4f * x² + %.4f * x + %.4f\n", a, b, c);

        // 计算剩余平方和和R²
        double sse = 0, sst = 0, mean_y = sum_y / n;
        for (int i = 0; i < n; i++) {
            double y_pred = a * x[i] * x[i] + b * x[i] + c;
            sse += pow(y[i] - y_pred, 2);
            sst += pow(y[i] - mean_y, 2);
        }
        double r_squared = 1.0 - (sse / sst);

        printf("Sum of squared errors: %.4f\n", sse);
        printf("Coefficient of determination (R²): %.4f\n", r_squared);
    }
    else if (t == 3) { // 三次函数拟合 y = ax³ + bx² + cx + d
        printf("Performing cubic fitting (y = ax³ + bx² + cx + d)...\n");

        double sum_x = 0, sum_y = 0, sum_xx = 0, sum_xy = 0;
        double sum_xxx = 0, sum_xxxx = 0, sum_xxy = 0;
        double sum_xxxxx = 0, sum_xxxxxx = 0, sum_xxxy = 0;

        for (int i = 0; i < n; i++) {
            double xi = x[i], yi = y[i];
            sum_x += xi;
            sum_y += yi;
            sum_xx += xi * xi;
            sum_xy += xi * yi;
            sum_xxx += xi * xi * xi;
            sum_xxxx += xi * xi * xi * xi;
            sum_xxy += xi * xi * yi;
            sum_xxxxx += xi * xi * xi * xi * xi;
            sum_xxxxxx += xi * xi * xi * xi * xi * xi;
            sum_xxxy += xi * xi * xi * yi;
        }

        // 构建线性方程组
        double A[4][5] = {
            {n,      sum_x,  sum_xx,  sum_xxx,    sum_y},
            {sum_x,  sum_xx, sum_xxx, sum_xxxx,   sum_xy},
            {sum_xx, sum_xxx, sum_xxxx, sum_xxxxx, sum_xxy},
            {sum_xxx, sum_xxxx, sum_xxxxx, sum_xxxxxx, sum_xxxy}
        };

        // 高斯消元法解方程组
        for (int k = 0; k < 4; k++) {
            double pivot = A[k][k];
            if (fabs(pivot) < 1e-10) {
                printf("Error: Singular matrix encountered\n");
                free(x);
                free(y);
                return;
            }

            for (int j = k; j < 5; j++) {
                A[k][j] /= pivot;
            }

            for (int i = 0; i < 4; i++) {
                if (i != k) {
                    double factor = A[i][k];
                    for (int j = k; j < 5; j++) {
                        A[i][j] -= factor * A[k][j];
                    }
                }
            }
        }

        double d = A[0][4];
        double c = A[1][4];
        double b = A[2][4];
        double a = A[3][4];

        printf("Fitted parameters: a = %.6f, b = %.6f, c = %.6f, d = %.6f\n", a, b, c, d);
        printf("Fitted equation: y = %.4f * x³ + %.4f * x² + %.4f * x + %.4f\n", a, b, c, d);

        // 计算剩余平方和和R²
        double sse = 0, sst = 0, mean_y = sum_y / n;
        for (int i = 0; i < n; i++) {
            double y_pred = a * x[i] * x[i] * x[i] + b * x[i] * x[i] + c * x[i] + d;
            sse += pow(y[i] - y_pred, 2);
            sst += pow(y[i] - mean_y, 2);
        }
        double r_squared = 1.0 - (sse / sst);

        printf("Sum of squared errors: %.4f\n", sse);
        printf("Coefficient of determination (R²): %.4f\n", r_squared);
    }
    else {
        printf("Error: Invalid function type t=%d (must be 1, 2, or 3)\n", t);
    }

    free(x);
    free(y);
    printf("=== Fitting completed ===\n");
}

// 可视化函数
void visualize(const char* filename) {
    printf("\n=== Visualizing dataset ===\n");

    // 构建完整路径
    char full_path[256];
    sprintf(full_path, "%s%s", DATA_DIR, filename);
    printf("Visualizing file: %s\n", full_path);

    FILE* file = fopen(full_path, "r");
    if (!file) {
        perror("Error opening file for visualization");
        return;
    }

    printf("Sample points (x, y):\n");
    double x, y;
    int count = 0;
    while (fscanf(file, "%lf %lf", &x, &y) != EOF && count < 10) { // 只显示前10个点避免刷屏
        printf("* (%.0f, %.2f)\n", x, y);
        count++;
    }
    if (count < 10) {
        printf("... (all %d points shown)\n", count);
    }
    else {
        printf("... (showing first 10 of %d points)\n", count);
    }
    fclose(file);
}

int main(int argc, char* argv[]) {
    printf("\n=== Curve Fitting Program ===\n");

    // 打印所有接收到的参数
    printf("Received %d arguments:\n", argc);
    for (int i = 0; i < argc; i++) {
        printf("argv[%d] = %s\n", i, argv[i]);
    }

    if (argc < 2) {
        printf("Usage:\n");
        printf("To generate data:\n");
        printf("  gendat <filename> <t(1-3)> <a> <b> [<c> <d>]\n");
        printf("    t=1: linear (y=ax+b) - requires a, b\n");
        printf("    t=2: quadratic (y=ax²+bx+c) - requires a, b, c\n");
        printf("    t=3: cubic (y=ax³+bx²+cx+d) - requires a, b, c, d\n");
        printf("\nTo fit data:\n");
        printf("  fitting <filename> ls <t(1-3)>\n");
        printf("    ls: least squares method\n");
        printf("    t: function type to fit (1, 2, or 3)\n");
        return 1;
    }

    if (strcmp(argv[1], "gendat") == 0) {
        if (argc < 5) {
            printf("Error: insufficient arguments for gendat\n");
            printf("Usage: gendat <filename> <t(1-3)> <a> <b> [<c> <d>]\n");
            return 1;
        }

        const char* filename = argv[2];
        int t = atoi(argv[3]);

        if (t < 1 || t > 3) {
            printf("Error: t must be 1, 2, or 3\n");
            return 1;
        }

        double a = atof(argv[4]), b = atof(argv[5]);
        double c = 0.0, d = 0.0;

        if (t >= 2 && argc < 6) {
            printf("Error: quadratic function requires c parameter\n");
            return 1;
        }
        if (t >= 2) c = atof(argv[6]);

        if (t == 3 && argc < 7) {
            printf("Error: cubic function requires d parameter\n");
            return 1;
        }
        if (t == 3) d = atof(argv[7]);

        printf("Generating dataset with parameters:\n");
        printf("  filename: %s\n", filename);
        printf("  function type: %d\n", t);
        printf("  a: %.2f, b: %.2f", a, b);
        if (t >= 2) printf(", c: %.2f", c);
        if (t == 3) printf(", d: %.2f", d);
        printf("\n");

        generate_dataset(filename, t, a, b, c, d);
        visualize(filename);
    }
    else if (strcmp(argv[1], "fitting") == 0) {
        if (argc < 4) {
            printf("Error: insufficient arguments for fitting\n");
            printf("Usage: fitting <filename> ls <t(1-3)>\n");
            return 1;
        }

        const char* filename = argv[2];
        const char* method = argv[3];
        int t = (argc >= 5) ? atoi(argv[4]) : 1;  // 默认t=1

        if (t < 1 || t > 3) {
            printf("Error: t must be 1, 2, or 3\n");
            return 1;
        }

        printf("Fitting parameters:\n");
        printf("  filename: %s\n", filename);
        printf("  method: %s\n", method);
        printf("  function type: %d\n", t);

        if (strcmp(method, "ls") == 0) {
            least_squares_fit(filename, t);
        }
        else {
            printf("Error: unknown fitting method '%s'\n", method);
            return 1;
        }
    }
    else {
        printf("Error: unknown command '%s'\n", argv[1]);
        return 1;
    }

    printf("\n=== Program completed ===\n");
    return 0;
}